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Abstract: We study identification of stochastic Wiener dynamic systems using so-called indirect 
inference. The main idea is to first fit an auxiliary model to the observed data and then in a 
second step, often by simulation, fit a more structured model to the estimated auxiliary model. 
This two-step procedure can be used when the direct maximum-likelihood estimate is difficult 
or intractable to compute. One such example is the identification of stochastic Wiener systems, 
i.e., linear dynamic systems with process noise where the output is measured using a non¬ 
linear sensor with additive measurement noise. It is in principle possible to evaluate the log- 
likelihood cost function using numerical integration, but the corresponding optimization problem 
can be quite intricate. This motivates studying consistent, but sub-optimal, identification 
methods for stochastic Wiener systems. We will consider indirect inference using the best linear 
approximation as an auxiliary model. We show that the key to obtain a reliable estimate is to use 
uncertainty weighting when fitting the stochastic Wiener model to the auxiliary model estimate. 
The main technical contribution of this paper is the corresponding asymptotic variance analysis. 
A numerical evaluation is presented based on a first-order finite impulse response system with 
a cubic non-linearity, for which certain illustrative analytic properties are derived. 

Keywords: system identification, indirect inference, identification of non-linear systems, 

Wiener systems, stochastic non-linear system. 


1. INTRODUCTION 

The idea of using two-stage identification methods, e.g., 
indirect inference, is by no means new in system identi¬ 
fication. Typically, a flexible auxiliary model is fitted to 
data, and then in a second step this estimated model is 
used to find a more structured model. A well-known ex¬ 
ample of such an approach is the indirect Prediction Error 
Minimization (PEM) method, Soderstrbm et al. (1991), 
where it is assumed that the model structure of interest 
can be embedded in a larger model structure for which 
the identification problem is more tractable. In a second 
step, the structured model is estimated from the larger 
model using a weighted non-linear least squares method. 
The indirect PEM method will, under certain assumptions, 
have the same asymptotic statistical properties as the 
Maximum-Likelihood (ML) method, but it can be more 
efficiently calculated. Indirect PEM is a special case of 
indirect inference, which was introduced in econometrics 
in Gourieroux et al. (1993). Their main motivation was 
identification problems for which the ML method is in¬ 
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tractable. They also proposed the use of Monte Carlo 
simulations to generate the cost to be minimized in the 
second step. 
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Fig. 1. Stochastic Wiener System 

The concept of indirect inference was introduced to the 
system identification community by Welsh et al. (2009) 
and Larsson et al. (2010). The aim of the current paper 
is to provide further insights into the potential use of 
indirect inference for the identification of scalar discrete 
time stochastic Wiener systems, illustrated in Figure 1, of 
the form 

z{t) = G(q)u(t) + v(t), 

y(t) = f(z(t)) + e(t), (1) 

with a stable transfer function G(q) (where q denotes the 
shift operator), an input signal sequence (n(f)}, white 
stationary process noise {i;(f)} with zero mean and vari¬ 
ance dy , an output signal sequence {y(t)}, and additive 
white stationary measurement noise (e(f)} with zero mean 
and variance a\. The input signal u(t) is assumed to 














be independent of the noises (the open-loop case), and 
may be a realization of a stochastic process with known 
probability distribution. To simplify the presentation, we 
will assume that the noise processes are independent and 
normal (gaussian) distributed. The more general case with 
coloured process noise can be handled using a predictor 
model. 

The main challenge is the non-linear function /(•), which 
means that we have a non-linear stochastic system where 
the process noise (i’(t)} propagates through a non-linear 
device. This typically corresponds to a non-linear sensor. 

1.1 Contribution 

The main objective of this paper is to study the indirect 
inference method using the best linear approximation as an 
auxiliary model and compare this to maximum-likelihood 
and prediction error minimization methods for identifica¬ 
tion of a stochastic Wiener system. The main challenge 
is to handle the non-linear process noise contribution. It 
is in principle possible to calculate the likelihood cost, or 
some approximation of it, using for example particle filters 
and/or Markov Chain Monte Carlo methods, but even for 
simple examples the maximum-likelihood approach leads 
to rather involved computations. A problem for prediction 
error minimization methods is that the probability density 
function of the conditional mean prediction error is quite 
complicated and involves convolution integrals. This is the 
reason for using more ad hoc identification techniques. A 
common approach is to use linear and gaussian approxi¬ 
mations. The main technical contributions are: 

• To connect the use of best linear approximation of 
stochastic Wiener systems with the method of indi¬ 
rect inference. We present the corresponding variance 
uncertainty weighting, which for this case includes the 
input signal. We also derive the asymptotic variance 
expression for the estimated model parameters. 

• Illustration of the results on a first-order finite im¬ 
pulse response model, for which it is possible to 
find analytic expressions. The statistical performance 
on this example of indirect inference based on best 
linear approximation is comparable to that that of 
maximum-likelihood estimation and prediction error 
minimization. The computations for the indirect in¬ 
ference method are just a fraction of the ones for 
calculating the maximum-likelihood estimate. The 
example also shows that the cost of using a non¬ 
linear sensor is increased uncertainty in the estimated 
model. 

The statistical theory for identification of stochastic 
Wiener systems is by no means complete, and the aim of 
this paper is to provide insights in some open important 
problems. 

1.2 System Identification 

Given measurements of the input and output signals 
{u(t),y(t)}, t = 1,..., N, the task is to identify a model 
of the stochastic Wiener system of the form, 
z(t) = G(q,9)u(t) +v(t), 

V(t ) = f(z(t)) + e{t) (2) 


where the model is parameterized by the parameter vector 
9 G R”. We assume that the true system can be described 
by 9 0 . The noise processes (n(t)} and (e(i)} are assumed to 
be independent normal distributed (gaussian) zero mean 
white noise. The corresponding noise variances <j% and 
<7g are assumed to be known, but could be added to 
the parameters to be estimated. This can, however, cause 
identifiability problems. We will study the case when the 
non-linear function /(•) is known. It is possible to extend 
our results to the case when the function /(•) also is 
estimated, which can result in an identifiability problem. 
The reason for these simplifications is to focus on the 
stochastic part due to the process noise v(t). 

Identification of Wiener systems is a well studied topic, see 
for example Greblicki (1992); Wigren (1993); Bai (2003); 
Zhu (2002); Enqvist and Ljung (2005); Pillonetto (2013); 
Giri and Bai (2010) and the references therein. It forms 
the basis for the identification of more general non-linear 
block diagram based models. However, many algorithms 
assume no process noise, which leads to a non-linear least 
squares problem minimizing the output error, i.e., the dif¬ 
ference between the measured and the simulated outputs. 
The maximum-likelihood method for stochastic Wiener 
systems was introduced in Hagenblad and Ljung (2000) 
and analysed in more detail in Hagenblad et al. (2008). 
The expectation-maximization algorithm for maximum- 
likelihood identification and the use of the particle filter 
have been studied in Wills et al. (2013); Wills and Ljung 
(2010). As recently pointed out in Wahlberg et al. (2014), 
the stochastic Wiener system identification problem can 
be viewed as a non-linear errors-in-variables problem, with 
the well-known bias problem due to input noise. This 
paper also proposes a prediction error minimization frame¬ 
work for identification of stochastic Wiener systems. The 
idea is to use the conditional mean predictor and notice 
that the variance of the prediction errors may be highly 
dependent on the input signal u(t). Hence, variance un¬ 
certainty weighting is most important in order to obtain 
reliable estimates. 

1.3 Outline 

ML and PEM methods for stochastic Wiener system iden¬ 
tification are reviewed in Sec. 2. In Sec. 3, we study how to 
use the indirect inference approach for the identification of 
stochastic Wiener models by using the BLA as the auxil¬ 
iary model. A first order FIR model example is outlined in 
Sec. 4, and the corresponding numerical study is presented 
in Sec. 5. The paper is concluded in Sec. 6. 

2. ML AND PEM 


As shown in Hagenblad and Ljung (2000), the negative log- 
likelihood function, given data and the normal distributed 
noise model (2), equals 

!• OO 

l(9) = -J2 l °S e- £{t ’ e ’ z) dz, 

+ _i J —OO 


t= 1 


where 


£(t,0,z) = 


lv(t) - /(A)] 2 , [z-G(q,e)u(t)f 


2 cr| 2cr2 

The ML estimate of 9 is obtained by minimizing 1(9). 
There are at least two challenges with the ML method for 



stochastic Wiener systems. First, to evaluate the negative 
log-likelihood cost at a certain parameter value 6 we have 
to calculate N integrals. This can be done rather efficiently 
using numerical integration and parallel computations. A 
difficulty is the integrand, where typically 

p -£(t,9,x) _ f U x small, 

~ 1 0, x otherwise. 

The example to be considered in Sec. 4 corresponds to 
£ (x) ~ x °, which means that the integrand decreases very 
rapidly to zero. 


The PEM approach avoids the exponential function in¬ 
tegration issue by using a weighted least squares cost 
function, see Wahlberg et al. (2014). The conditional mean 
predictor of y(t) for given u(t) and 9 is 

y(t, 0) = E v {f(G(q, 0)u{t) + v(t))}. (3) 

Notice that the prediction error variance depends on the 
input signal u(t). The optimally weighted quadratic PEM 
cost-function, see Wahlberg et al. (2014), is 


Vn(9) 


J_ g2 (C 9 ) 


(4) 


with prediction error e(t,9) = y(t) — y(t,9). The variance 
weighting is calculated at a consistent initial estimate 
Or, e.g., the PEM estimate without weighting. The use 
of weighting is important to obtain reliable estimates 
and depends here on the input signal u(t). The PEM 
estimate based on (4) is not asymptotically efficient for 
non-linear functions, since we use a weighted quadratic 
cost-function. However, by using a cost-function based on 
the probability density function of e(t,9), we obtain an 
asymptotically efficient PEM estimate, see Ljung (1999). 
The computations will then be similar to the ML case, 
involving multiple integral calculations with exponential 
functions. 


3. INDIRECT INFERENCE USING BEST LINEAR 
APPROXIMATION 

3.1 Best Linear Approximation of Stochastic Wiener 
Systems 

Let us illustrate the concept of indirect inference applied 
to a stochastic Wiener system with a known non-linear 
function. A common first approach is to fit a linear model 
to the observed data. It is well known that if the input 
signal is normal (gaussian) distributed, then the Best 
Linear Approximation (BLA), see Ljung (2001) or (Giri 
and Bai, 2010, Chapter 13), is a scaled version of the linear 
dynamics transfer function G(q) of the Wiener system. 
It is perhaps less well known that the same result holds 
for stochastic Wiener systems with gaussian process noise 
v(t). This extension follows more or less from Bussgang’s 
theorem, Bussgang (1952): If z(t) is a normal distributed 
stationary stochastic process with zero mean and if the 
non-linear transformed process y(t) = f(z(t)) has zero 
mean, then 

E {y(t)z(t - r)} = b Q E{z(t)z(t - r)}, b 0 = E{f'(z(t))}. 

Adding independent gaussian process noise to the input 
signal contribution, z(t) = G(q)u(t) +v(t), makes z(t) still 
normal distributed and Bussgang’s theorem holds. Now we 


are only interested in computing the partial correlations 
using the relations 

E {y(t)u(t - t)} = b 0 E{z(t)u(t - t)} 

= b 0 E{[G(q)u(t)]u(t - t)}. (5) 

A recent proof of (5) can be found in Banelli (2012), but 
it already follows from Nutall (1958). This result proves 
that the BLA of a stochastic Wiener system with a normal 
distributed input signal again equals 

GsLA{q) = boG(q) 

i.e., a scaled version of the linear transfer function. Notice 
that b 0 now also depends on the statistics of the process 
noise v as well as of the input u. This result can be gen¬ 
eralized to separable stochastic processes, Nutall (1958); 
Enqvist and Ljung (2005). 

The BLA of a stochastic Wiener system for more general 
input signal sequences can be obtained by simulations or 
by analytic calculations (if the distribution of the input 
signal is known). 

3.2 An Optimally Weighted Indirect Inference Algorithm 

We will now describe how to apply indirect inference to 
the stochastic Wiener system identification problem using 
a two-step procedure based on BLA. From the given data 
{u(t),y(t)}, t = 1,..., N, form the BLA cost-function 

1 N 

Qn( 0) = - Gun(q , P)u(t)} 2 . (6) 

Here we have used the notation Gu n (q, 0) to stress that 
this is a linear transfer function parameterized by (3 £ R m . 
It is in general of higher order compared to G(q , 9) in (2), 
that is to is typically larger than n, the dimension of 9. 
We have used an Output Error PEM cost-function (6), 
but it is also possible to use a more general PEM model 
structure, as described in Ljung (2001). 

Step 1: Identify the BLA within the model structure 
Gun(q,f3), 0 6 R m , of the system from given data using 
the cost-function (6) 

0n = arg irdn{Qjv(/3)}, 

i.e., the BLA estimate equals GsLAiq) — Gu n {q, 0 n)- 

The next question is to figure out the functional relation 
0(9), R m i—^ R" that describes how the auxiliary estimate 
0n asymptotically depends on the underlying structured 
model of interest, as a function of the model parameter 
vector 0. Now assume that Qn(0) converges (almost 
surely) as the number of data N —> oo to Q(0, 9 0 ), which 
depends on the true system parameter vector 0 o . Here 
Q(0,9)=E{(y(t)-G lin (q,0)u(t)] 2 } 

= E v , u {[f(G(q, 9)u(t) + v(t)) - Gii„ (q, 0)u(t)] 2 }. 
The expectation is with respect to both the input signal 
and the process noise. It is also be possible to only use 
expectation with respect to the process noise for a given 
input sequence. Let 

0(9) = arg min {Q(0,9)}, 

and use the notation 

0o = 0(9 o), 

for the corresponding true parameter vector. This leads to: 



Step 2: (Analytic) Estimate the structured model pa¬ 
rameter vector 0 by solving 

On = argmin[/?(#) - p N ] T W[P{0) - p N ], (7) 

0 

where IE is a positive definite weighting matrix to be 
specified. 

For certain examples of non-linear functions and distribu¬ 
tions it may not possible to analytically find the function 
/ 3(0 ). One can then resort to Monte Carlo simulations. Let 

Qn,s(P,0) = 
s N 

C E N *Mf) + Veit)) Gliniq , P)u(t)f}. 

D s=l JV t=1 

Here (i> s (f)}, t = 1 s = 1_ ,S, is a generated 

noise realization of the process noise v(t') and S is the total 
number of realizations used in the Monte Carlo Simulation. 
Let 

Pn,s(0) =argmm {Qn,s(P,0)}. 

Step 2: (Simulated) Estimate the structured model 
parameter vector 0 by solving 

On,s = argmin[/3 JV> s(6») - Pn] T W[Pn,s(0) - Pn], (8) 

u 

where W is a positive definite weighting matrix to be 
specified. 

The optimal weighting matrix, W, equals the inverse of the 
covariance matrix of the auxiliary estimate Pn, compare 
Expression (9.11) for the asymptotic covariance matrix in 
Ljung (1999). In our setting this translates to 

W = (9) 

/ 0 = lim Co v{VN^^-(P 0 ,d 0 )}, Jo = ^M(Po,0 o ). 

N^oo Op Op z 

In practise, one has to use a consistent estimate of W. 

3.3 Performance Analysis 

Next, we determine the asymptotic properties of the final 
estimate 9n- This can also be done using Taylor series 
expansion arguments as described in Complement C4.4 
in Soderstrom and Stoica (1989). The properties of the 
function P(0), R m R ra , are very important in order to 
obtain a consistent estimate. It should be possible to invert 
this function, 6 = a(P(9)), i.e., a(-) is a left inverse of P(-). 
Define the Jacobian matrix 

G=|j(d 0 )e R mxn , (10) 

and recall that we use the weighting matrix 

W = [ lim Cov{VN[P n - Po}}}^ 1 . (11) 

N—> oo 

Variance Expression: The asymptotic covariance matrix of 
On defined by (7) equals 

lim Co v{Vn[§ n - Oo]} = [G t WG]-\ (12) 

N—>oo 

where G is defined by (10) and W by (11). 

To prove this result, we will make a Taylor series expansion 
of the cost-function at the minimizing 0: 

V w (0) = ]P(0) - p N ] T W[P(0) - Pn] =► V^0 N ) = 0, 

which gives [On — 0 o ] ~ —\Yw{p o )]~ 1 Vf l ,{0o). The Hessian 
Vff(0 o ) has to be invertible (positive definite for unique 


local minimum) for this to hold. As shown in Complement 

4.4 in Soderstrom and Stoica (1989) 

Vf,(0 o ) = 2 G t W[P n - p 0 ] + 0{1/N) 

V"(0 O ) = 2G t WG + <D(l/y/N). 

Hence 

[On - 0 o ] w -[G t WG]~ 1 G t W[Pn - Po]- 
Since limjv-j-oo Cov{\/iV[j8jv— Poi\ = IE -1 , this proves (12). 

When using the Monte Carlo simulation approach in 
Step 2, (8), the asymptotic covariance matrix of 0n,s 
should be amplified by the factor (1 + l/S 1 ), due to the 
extra uncertainty from the simulations, see Heggland and 
Frigessi (2004). 

3-4 Comments 

The suggestion to use the BLA as an auxiliary model is 
rather ad hoc, but a very common choice in recent meth¬ 
ods for identification of non-linear systems, Schoukens 
et al. (2005); Pintelon and Schoukens (2012); Schoukens 
et al. (2014); Schoukens and Rolain (2012); Sjoberg and 
Schoukens (2012). Ljung (2001) contains an overview of 
the role of BLA in system identification. Another option 
would be to use the biased estimate from minimizing 

N 

Z (^) = - f(G{q,p)u(t)] 2 

t=l 

and then use Step 2 to remove the bias of this estimate. 
A challenge is to find even more efficient auxiliary models. 
The key property is that p(0) should have a left inverse 
(identifiability) and 

G = ^(0 o )e R mxn (13) 

should be “large” (sensitive) and at the same time p should 
be easy to estimate from the given data. 

3.5 Indirect Inference 

As mentioned in the introduction the indirect inference 
approach was developed in Gourieroux et al. (1993) as 
a way to find consistent model estimates even when the 
ML method is intractable. Our description of the indirect 
inference approach applied to BLA in the previous section 
is mainly based on Heggland and Frigessi (2004). The 
theory of indirect inference is more general, and the key 
step is to choose the auxiliary model parameterized by P 
and the data-driven cost-function Qn{P)- A convergence 
analysis of the general indirect inference method presented 
in Heggland and Frigessi (2004) is based on using the 
concept of finite dimension auxiliary statistics. 

If Qn(P) is based on a sufficient statistics the indirect 
inference method is efficient (except for the factor (1 + 
1 //S')). This is not the case for stochastic Wiener systems, 
for which no finite dimensional sufficient statistics exists. 
For the BLA approach we are only using a second order 
statistics. 

4. ILLUSTRATING EXAMPLE 

We will use the following simple example to try to further 
understand the properties of the stochastic Wiener system 
identification methods described in the previous sections, 



z[t) = 9u(t) + u(t — 1) + v(t) 
y(t) = [z(t)] 3 + e(t). (14) 

The same example was used in Larsson et al. (2010); 
Wahlberg et al. (2014). We will study two different input 
white noise distributions, a normal (gaussian) and an 
uniform distribution, both with variance cr 2 . The noises 
are assumed to be white zero mean normal distributed 
with variances cr 2 and cr 2 , respectively. 


The maximum-likelihood cost-function 


N 


m 


= ~H logE^Tff 

t=l *2^ 


a v --k I [y(t)-(9u(t)+u(t-l)+a v v) 3 ]' i 

e 2<T e } 


is calculated using a Gauss Hermite approximation of 
order 1000. The reason for this high order is the e~ x tends 
to zero very quickly and the integral is rather difficult to 
approximate. We have noticed that increased process noise 
variance a 2 makes it even more challenging. 


The conditional mean predictor equals 

y(t, 9) = 9 3 u 3 {t) + u 3 (t — 1) + 3 (9u 2 (t)u(t — 1) 

+ 9u(t)u 2 {t - 1) + 9u{t)cr 2 + u(t - l)cr 2 ), (15) 

and is used in the PEM method (4). 


In order to evaluate the indirect inference method we 
will study the BLA. Assume first that the input signal is 
normal distributed with zero mean and variance cr 2 . The 
BLA of model order one, i.e., based on minimizing 

E{(2/(^) - + #2 u(t - l)) 2 } 

is 

(|$)= [ 3<^ 2 + 3 « + 0 | (()'< ( 16 ) 

for which 



since the BLA for gaussian input signal is just a scaled 
version of G(q,9). This also gives that (3(a) has a left 
inverse and we have identifiability of 9 from /3. 


The corresponding BLA when the input signal is uniformly 
distributed with zero mean variance a 2 is 


f 0 i(g) \ _ ( + 3 * it “ + ) 

U ( 9 ) J- (^ + 3 + 

For the indirect inference method we will use 
1 N 

Qn{P) = jy - A«(t) + PM* -1)) 2 


(17) 


and hence 

9Qn = _2_ ^ ( ( y(t) - /3i u(t) + p 2 u(t - 1 ))u{t) \ 

dp N ^ V ( y{t ) - /3i«(i) + fi 2 u(t - 1 ))u(t - 1) ) 


and we need to calculate the weighting matrix W from (9). 
Here J 0 = cr 2 1 and the tedious work is to calculate I 0 . It 
can, however, be estimated as 

1 N 

= 77 [fezW _ ^ u (*) _ PM* - i)] 2 x 

V t= 1 

( U 2 {t) u(t)u(t — 1)7 

\u(t))u(t— 1) u 2 {t— 1) J 


Comment: The key to get the order one indirect inference 
method to work is to use the weighting W = I~ 1 . The 
paper Larsson et al. (2010) does not use weighting, which 
explains their non-intuitive simulation result that it is 
better to use a zero order BLA model than a first order 
one (which should contain more information). 

We will finally study the zero order model case, for which 
no weighting is needed. Here 

1 N 

Q n (P) = j-f ^2{y(t) - Piu{t)) 2 

v t=i 

with the BLA given by Pi(9) in (16) and (17). An alter¬ 
native cost-function is 

1 N 

Qn{/ 3) = ^2(y(t) - p 2 u{t - 2)) 2 , 

v t=1 

with the BLA given by p 2 (9) from (16) or (17). A problem 
here is that fi 2 ( 9 ) is quadratic in 9 and we do not have an 
unique solution with respect to 9. 

5. SIMULATION RESULT 

We will use the following numerical values for the system 
parameters in (14): 

9 0 = 0.5, (j 2 = 0.2, (j 2 = 0.1. 

We will use the analytic Step 2, (7), for the indirect 
inference based on BLA, and will evaluate the following 
methods for N = 1000 observations: 


Method 1 
Method 2 
Method 3 
Method 4 
Method 5 


ML 

PEM with optimal weighting 

Zero order indirect inference 

First order indirect inference, no weighting 

First order indirect inference, with weighting 


The following table summarizes the simulation results for 
the normal distributed input sequence with zero mean and 
variance a 2 = 1/3 over 1000 noise and input realizations: 


Method 1: mean: 
Method 2: mean: 
Method 3: mean: 
Method 4: mean: 
Method 5: mean: 


0.5025 

std: 0.0303 

0.5014 

std: 0.0349 

0.4983 

std: 0.0446 

0.4977 

std: 0.0554 

0.4982 

std: 0.0418 


The next table summarizes the simulation results for the 
uniform distributed input sequence with zero mean and 
variance <r 2 = 1/3 over 1000 noise and input realizations 


Method 1: mean: 
Method 2: mean: 
Method 3: mean: 
Method 4: mean: 
Method 5: mean: 


0.4994 

std: 0.0325 

0.4984 

std: 0.0346 

0.4988 

std: 0.0454 

0.4984 

std: 0.0458 

0.4987 

std: 0.0377 


The conclusion from the numerical study is that the ML, as 
expected, gives the best performance. However, PEM and 
indirect inference with optimal weighting give also good 






results. The computational times for these methods are 
only a fraction of that of ML. The simulation study also 
shows that a zero order model can give better results than 
using a first order model and no weighting, c.f., Larsson 
et al. (2010). 

It should be noted that the asymptotic variance when 
f(x) = x equals (crj + a),)/(afN), with a standard 
deviation equal to 0.0173 for the example above. Hence, 
the non-linear function f(x) = x 3 gives a more difficult 
estimation problem than the standard case f(x) = x. 

6. CONCLUSION 

In this paper, we have utilised the indirect inference 
method based on the best linear approximation to identify 
stochastic Wiener systems. The results show that to obtain 
a reliable estimate, it is important to use an optimal 
weighting when estimating the structured model from the 
auxiliary model. The weighting here is the inverse of the 
covariance matrix of the BLA parameter estimate. We 
have analyzed the statistical properties of the correspond¬ 
ing indirect inference estimate using results from system 
identification. The methods have been evaluated using a 
first order FIR system with a cubic non-linearity. The 
simulation results demonstrate that the proposed indirect 
inference BLA approach performs quite well compared to 
the ML and weighted PEM methods. A major advantage 
is that the indirect inference algorithms are very compu¬ 
tationally fast and direct to implement. 

There are many open questions when it comes to identi¬ 
fication of stochastic non-linear dynamic systems. For ex¬ 
ample, performance results to guide the choice of sensors. 
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